:- set(i,2).
:- set(nodes,20000).
:- set(noise,5).
:- set(c,3).
:- set(minpos, 2).

:- set(test_pos, './mutagenesis/188_friendly_aleph_lazy/t5.f').
:- set(test_neg, './mutagenesis/188_friendly_aleph_lazy/t5.n').

:- modeh(1,active(+drug)).


:- modeb(1,lumo(+drug,-number)).
    :- modeb(1,lumo(+drug,#number)).


:- modeb(1,logp(+drug,-number)).
    :- modeb(1,logp(+drug,#number)).

:- modeb(*,bond(+drug,-atomid,-atomid,#number)).
:- modeb(*,bond(+drug,+atomid,-atomid,#number)).
    :- modeb(*,bond(+drug,+atomid,-atomid,-number)).

:- modeb(*,atm(+drug,-atomid,#element,#number,-charge)).
    :- modeb(*,atm(+drug,-atomid,#element,#number,#charge)).
	:- modeb(*,atm(+drug,-atomid,-element,-number,#charge)).
	:- modeb(*,atm(+drug,-atomid,-element,-number,-charge)).


:- modeb(*,benzene(+drug,-ring)).
:- modeb(*,carbon_5_aromatic_ring(+drug,-ring)).
:- modeb(*,carbon_6_ring(+drug,-ring)).
:- modeb(*,hetero_aromatic_6_ring(+drug,-ring)).
:- modeb(*,hetero_aromatic_5_ring(+drug,-ring)).
:- modeb(*,ring_size_6(+drug,-ring)).
:- modeb(*,ring_size_5(+drug,-ring)).
:- modeb(*,nitro(+drug,-ring)).
:- modeb(*,methyl(+drug,-ring)).
:- modeb(*,anthracene(+drug,-ringlist)).
:- modeb(*,phenanthrene(+drug,-ringlist)).
:- modeb(*,ball3(+drug,-ringlist)).


:- modeb(*,member(-ring,+ringlist)).
:- modeb(1,member(+ring,+ringlist)).
:- modeb(1,connected(+ring,+ring)).



:- determination(active/1,atm/5).
:- determination(active/1,bond/4).

:- determination(active/1,lumo/2).
:- determination(active/1,logp/2).

:- determination(active/1,benzene/2).
:- determination(active/1,carbon_5_aromatic_ring/2).
:- determination(active/1,carbon_6_ring/2).
:- determination(active/1,hetero_aromatic_6_ring/2).
:- determination(active/1,hetero_aromatic_5_ring/2).
:- determination(active/1,ring_size_6/2).
:- determination(active/1,ring_size_5/2).
:- determination(active/1,nitro/2).
:- determination(active/1,methyl/2).
:- determination(active/1,anthracene/2).
:- determination(active/1,phenanthrene/2).
:- determination(active/1,ball3/2).
:- determination(active/1,member/2).
:- determination(active/1,connected/2).


:- modeb(1,gteq(+number,#number)).
:- modeb(1,lteq(+number,#number)).
:- modeb(1,geq(+number,#number)).
:- modeb(1,leq(+number,#number)).

:- lazy_evaluate(gteq/2).
:- lazy_evaluate(lteq/2).
:- lazy_evaluate(geq/2).
:- lazy_evaluate(leq/2).

:- determination(active/1,gteq/2).
:- determination(active/1,lteq/2).
:- determination(active/1,geq/2).
:- determination(active/1,leq/2).

% definition to use during lazy evaluation
leq([P,N],Value):-
        !,
        find_cutoff(lte,P,N,Value).
% definition to use during normal evaluation
leq(X,Value):-
	number(X), number(Value), !,
        X < Value.
% definition to use during construction of bottom clause
leq(X,X).


% definition to use during lazy evaluation
geq([P,N],Value):-
        !,
        find_cutoff(gte,P,N,Value).
% definition to use during normal evaluation
geq(X,Value):-
	number(X), number(Value), !,
        X > Value.
% definition to use during construction of bottom clause
geq(X,X).


% definition to use during lazy evaluation
lteq([P,N],Value):-
        !,
        find_cutoff(lte,P,N,Value).
% definition to use during normal evaluation
lteq(X,Value):-
	number(X), number(Value), !,
        X =< Value.
% definition to use during construction of bottom clause
lteq(X,X).


% definition to use during lazy evaluation
gteq([P,N],Value):-
        !,
        find_cutoff(gte,P,N,Value).
% definition to use during normal evaluation
gteq(X,Value):-
	number(X), number(Value), !,
        X >= Value.
% definition to use during construction of bottom clause
gteq(X,X).



find_cutoff(gte,[],N,Value):-
	!,
	list_max(N,Value).
find_cutoff(gte,P,[],Value):-
	!,
	list_min(P,Value).
find_cutoff(lte,[],N,Value):-
	!,
	list_min(N,Value).
find_cutoff(lte,P,[],Value):-
	!,
	list_max(P,Value).
find_cutoff(C,P,N,Value):-
        sort(P,Ps), sort(N,Ns),         	% Ps, Ns are in ascending order
        length(P,P0), length(N,N0),
        PriorFrac is P0/(P0+N0),
        PriorBits is -log(PriorFrac)/log(2),    % bits required to encode a prior pos
        find_loss(C,Ps,P,N,P0,N0,PriorBits,PLosses),   % losses made by each pos cutoff
        find_loss(C,Ns,P,N,P0,N0,PriorBits,NLosses),   % losses made by each neg cutoff
	keysort(PLosses,[PL1-PV1|_]),
	keysort(NLosses,[NL1-NV1|_]),
	(NL1 < PL1 -> Value is NV1 ; Value is PV1).


find_loss(_,[],_,_,_,_,_,[]).
find_loss(C,[X|T],Ps,Ns,P0,N0,PriorBits,[Loss-X|T1]):-
	find_covered(C,X,Ps,P0,P1),
	find_covered(C,X,Ns,N0,N1),
        PostFrac is (P1+0.001)/(P1+N1),
        PostBits is -log(PostFrac)/log(2),      % bits required to encode a post pos
        GainBits is P1*(PriorBits - PostBits),  % total bits gained
        Loss is -GainBits,
        find_loss(C,T,Ps,Ns,P0,N0,PriorBits,T1).

find_covered(lte,_,[],N,N):- !.
find_covered(lte,X,[Y|T],N,N1):-
	Y > X, !,
	length([Y|T],N0),
	N1 is N - N0.
find_covered(gte,_,[],_,0):- !.
find_covered(gte,X,[Y|T],_,N1):-
	Y >= X, !,
	length([Y|T],N1).
find_covered(C,X,[_|T],N,N1):-
	find_covered(C,X,T,N,N1).

list_max([H|T],X):-
	find_elem(max,T,H,X).

list_min([H|T],X):-
	find_elem(min,T,H,X).

find_elem(_,[],X,X).
find_elem(max,[H|T],X,Y):-
	(H > X -> find_elem(max,T,H,Y); find_elem(max,T,X,Y)).
find_elem(min,[H|T],X,Y):-
	(H < X -> find_elem(min,T,H,Y); find_elem(min,T,X,Y)).	



% type information

drug(D):-
	name(D,[_|X]), name(Num,X), int(Num),
	Num >= 1, Num =< 230, !.

atomid(A):-
	name(A,[_|X]),
	append(Z,[95|Y],X),
	name(N1,Y),
	name(N2,Z),
	int(N1), int(N2),
	N2 >= 1, N2 =< 230,
	N1 =< 500, !.

append([],A,A).
append([H|T],A,[H|T1]):-
        append(T,A,T1).



charge(X):-
        float(X).
energy(X):-
        float(X).
hydrophob(X):-
        float(X).

% ring([_|_]).

ring(X):-
	name(X,[114|_]).

ringlist(X):- 
	name(X,[114|_]).

% ringlist([]).
% ringlist([Ring|Rings]):-
        % ring(Ring),
        % ringlist(Ringlist).



element(br).
element(c).
element(cl).
element(f).
element(h).
element(i).
element(n).
element(o).
element(s).


% background knowledge 

gteq(X,Y):-
	not(var(X)), not(var(Y)),
	float(X), float(Y), 
	X >= Y, !.
gteq(X,X):-
	not(var(X)),
	float(X).

lteq(X,Y):-
	not(var(X)), not(var(Y)),
	float(X), float(Y),
	X =< Y, !.
lteq(X,X):-
	not(var(X)), 
	float(X).




:- consult(atom_bond).
:- consult(ind1).
:- consult(inda).
:- consult(log_mutag).
:- consult(logp).
:- consult(lumo).
:- consult(ring_struc).
:- consult(ring_theory).

